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ABSTRACT 

We present a new code for performing general-relativistic radiation-hydrodynamics simula- 
tions of accretion flows onto black holes. The radiation field is treated in the optically-thick ap- 
proximation, with the opacity contributed by Thomson scattering and thermal bremsstrahlung. 
Our analysis is concentrated on a detailed numerical investigation of hot two-dimensional, 
Bondi-Hoyle accretion flows with various Mach numbers. We find significant differences 
with respect to purely hydrodynamical evolutions. In particular, once the system relaxes to 
a radiation-pressure dominated regime, the accretion rates become about two orders of mag- 
nitude smaller than in the purely hydrodynamical case, remaining however super-Eddington 
as are the luminosities. Furthermore, when increasing the Mach number of the inflowing gas, 
the accretion rates become smaller because of the smaller cross section of the black hole, but 
the luminosities increase as a result a stronger emission in the shocked regions. Overall, our 
approach provides the first self-consistent calculation of the Bondi-Hoyle luminosity, most of 
which is emitted within r ~ 100M from the black hole, with typical values L/L^dd — 1 — 7, 
and corresponding energy efficiencies r] Bn ~ 0.09 — 0.5. The possibility of computing lumi- 
nosities self-consistently has also allowed us to compare with the bremsstrahlung luminosity 
often used in modelling the electromagnetic counterparts to supermassive black-hole binaries, 
to find that in the optically-thick regime these more crude estimates are about 20 times larger 
than our radiation-hydrodynamics results. 
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1 INTRODUCTION 

Numerical relativity faces an embarrassing gap between the ac- 
curacy with which it computes the gravitational-wave emission 
from the dyn amics of compact objects such as black holes and 
neutron stars dMcWilliamsll2010l: lDuezll2010l : iRezzolla et a"i1l201ll : 



ISekiguchi et al. 201 lh and the very rough estimates of the elec- 



tromagnetic emission that can be currently com 


puted with state 


of the art numerical codes (Farrisetal. 


20081: 


Palenzuela et al. 


l2009l:lM6sta et alJ20ld:lBode et alJ2O10t 


Zanotti et al.ll2010h. The 



strongest limitation preventing a more realistic description of the 
emitted electromagnetic radiation is the modelling of radiative 
transfer in the gas, which is often neglected in relativistic calcula- 
tions in view of the large computational costs involved. This prob- 
lem is of course common to a large class of relativistic simulations, 
but it becomes particularly apparent in those cases where an accu- 
rate computation of the emitted luminosity is at least as important 
as providing a faithful description of the dynamics. Among such 
cases, accretion onto compact objects is perhaps the most impor- 
tant one. 



Until a few years ago, the time dependent solution of the rel- 
ativistic radiation hydrodynamics equations of accretion flows was 
performed in one spatial dimension only and, typically, through 
Lagrangian finite-difference schemes or through the so called lin- 
earized bl ock-implicit algori t hms. Starti ng from the p ioneering 
works by Gilden & Wheeler (1980) and IVitellol dl984h . relevant 
results concerning spherical accretio n onto black holes we re ob- 
tained with a Lagrangian code by IZampieri et al.l dl996L who 
were able to solve the radiation hydrodynamics equations both 
in the optically thin and in the optically thick regime, by means 
of the projected s ymmetr i c trac e-free (PSTF) moment formal - 
is m introduced by |Th orne dl98lh and subsequently reformulated 
by IRezzolla & Milled i 1994 ) for spherical flows. Such formalism 
provides one of the most accurate approximations to the solution 
of the radiation transfer equations, and, in analogy to what is done 
in fluid dynamics, it allows to define moments of the radiation field 
similarly to how density, momentum and pressure of a medium are 
defined as moments of the distribution function. As a result, instead 
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of following rays, the moment equations are solved directly with an 
Eulerian or a Lagrangian codejj 

Despite these initial efforts, the time-dependent solution of 
the relativistic radiation-hydrodynamics equations in more than one 
spatial dimensions remains very challenging. Nowadays, the mul- 
tidimensional numerical codes available can be divided in two ma- 
jor classes, accounting for, separately, the optically thin regime 
or the optically thick one. The former class is mainly focused 
on providing a realistic modelling of core-collapse supernovae, 
by employing Boltzmann neutrino transport, state-of- the- art neu- 
trino interactions, and general relativi ty. Relevant achievement s 
have been obtained over the years by [ .Mezzacap pa et al.l (120 1 ) ; 
iBruenn et aD(l200lh:llTebendorfer et al.ld200ll) : lLiebendorfer et all 
d2005h ;lMesser et al who, among other things, showed the 

importance of multidimensional simulations to model the shock re- 
vival via neutrinos in a supernova explosion. In a different physical 
context, namely that of accr etion discs around b lack holes, but still 
in the optically thin regime, iNoble et all d2009h considered an ap- 
proximate treatment in which radiation is described through a loss 
term in the energy equation. They used fully relativistic ray-tracing 
techniques to compute the luminosity received by distant observers. 
For a disc with aspect ratio H/r ~ 0.1 accreting onto a black hole 
with spin parameter a = 0.9, they found a si gnificant dissipation 
beyon d that predicted by the classical model by Novikov & Thornd 
dl973h . 

The numerical investigation of the optically thick regime, on 
the other hand, has rec eived less attention. The seminal work by 
iHsieh & Spiegel dl976h already set the basis for the formulation 
of the relativistic radiation-hydrodynamics equations in conserva- 
tion form and therefore suitable for an Eulerian numerica l imple 



menta ti on. La t er on, interestin g advan c es we re obtained by Shapiro 



19961) . IParkl d2006h and iTakahashil d2Q07h . Finally, iFarris et al 



(2008) have shown that for optically-thick gases and gray-body 
opacities, the general relativistic radiation-hydrodynamics equa- 
tions can indeed be written in conservation form, thus allowing 
for the use of numerical methods based on Riemann solvers that 
have been successfully adopted by many relativistic hydrodynam- 
ics and magnetohydrodynam ics codes. Very recen tly, and while this 
paper was being completed, IShibata et all (EoT 1) have presented a 
modified truncated moment formalism allowing for the conserva- 
tive formulation of the relativistic radiation-hydrodynamics equa- 
tions both in the optically thin and in the optically thick limit. This 
formulation could represent a major step forward with respect to 
prese nt leakage schemes accounting for the f ree streaming of radi- 
ation (Sekiguchi 2010; Se kiguchi et all201lh . 

I n this paper, which is th e first of a series, we extend our ECHO 
code Zann a et al ] l2007h by following the strategy suggested 
by IFarris et ah I d2008l) . and concentrate on one of the simplest 
accretion flows scenarios, namely: Bondi-Hoyle accr etion onto a 
black hole. This problem has recently been studied by IFarris et all 
boioh in the context of merging supermassive black hole binaries 
in full general relativity, but neglecting the back-reaction of radia- 
tion onto matter. By assuming that opacity is made of contributions 
by Thomson scattering and thermal bremsstrahlung, we compute 
here the luminosity emitted in hot Bondi-Hoyle accretion onto a 



1 In a no nrelativistic context, recent i nter esting developments have be en re- 
ported bv lPetkova & Springell d20Q9l) and lPetkova & Springell feOld) . who 
adopted the moment formalism within the SPH Gadget code and used a 
variable Eddington tensor as a closure rel ation, following t he Op tically Thin 
Variable Eddington Tensor suggestion of Gnedi n & Ab el ( 200lh . 



black hole. As the flow relaxes to a radiation-pressure dominated 
regime, we find significant differences with respect to purely hy- 
drodynamical evolutions. In particular, the accretion rates drop of 
about two orders of magnitude when compared to the purely hy- 
drodynamical case, remaining however super-Eddington. Further- 
more, we find that larger inflow velocities lead to smaller accretion 
rates (because of the smaller cross section of the black hole) but to 
larger luminosities (because of the stronger emission in the shocked 
regions). 

The plan of the paper is the following: We first describe the 
numerical methods in Section [2] while the validation of the code is 
presented in Section [3] Section |4] on the other hand, is devoted to 
radiative Bondi-Hoyle accretion flows and contains the main results 
of our work. The conclusions are presented in Sec. [5] We assume 
a signature ( — ,+,+,+) for the spacetime metric and we will use 
Greek letters (running from to 3) for four-dimensional spacetime 
tensor components, while Latin letters (running from 1 to 3) will be 
employed for three-dimensional spatial tensor components. More- 
over, we set c = 1, G = 10 -10 and extend the geometric units 
by setting m p /kB = 1, where m p is the mass of the proton, while 
ks is the Boltzmann constant. We have maintained c, G, and ks in 
a explicit form in those expressions of particular physical interest. 
Appendix [A] describes the extended geometrized system of units 
adopted in the code. 



2 RELATIVISTIC RADIATION HYDRODYNAMICS 
2.1 Covariant formulation 

The total momentum-energy tensor T a/3 of a fluid immersed in a 
radiation field comprises two terms: T a(3 = T^f + . The first 
one is the ordinary one describing the energy and momentum of the 
matter 



n a/3 



; a 8 i olB 

phu u + pg j 



(i) 



where g a ^ is the spacetime metric tensor, u a is the four-velocity 
of the fluid, while p, h = 1 + e + p/p, e and p are the rest-mass 
density, the specific enthalpy, the specific internal energy, and the 
thermal pressure, respectively. All of these quantities are measured 
in the comoving frame of the fluid. The thermal pressure is related 
to p and e through an equation of state (EOS), and we will here 
consider an ideal-gas, for which the EOS is expressed as 



p = pe ( 7 - 1) , 



(2) 



where 7 is the (constant) adiabatic index of the gas. The sec- 
ond i _j£nri_des£ri^ field and is given 
by dMihalas & Mihaladl 1984 IShapiroll 19961) 



I„N a N p dvdn, 



(3) 



where I u = lv (x 01 , N l , v) is the specific intensit)Q of the ra- 
diation, N a is the four-vector defining the photon propagation 
direction, dv is the infinitesimal frequency and dQ is the in- 
finitesimal solid angle around the direction of propagation. We re- 
call that the direction of propagation of the photon is defined as 
N a = p a /hpiv, where p a is the photon four-momentum, while 
/ipi and v are, respectively, the Planck constant and the pho- 
ton frequency as measured in the comoving frame of the fluid. 



We note that I u is an energy flux per unit time, frequency and solid angle, 



so that in cgs units it has dimensions of erg cm 



Hz" 
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Since the two terms vdvdfl an d I^/is 3 are relativistic invari- 
ants fevbicki & Lightmanl [l986h . their product with the tensor 
p ol p /B is still a tensor, and indeed it provides the integrand of Eq. 

In the frame comoving with the fluid, the moments of the radi- 
ation field are the energy density, the radiation flux and the radiation 
stress tensor, which are respectively given by 



- / IjydvdQ, 
c J 

h ">! 
-II 



I u disdQN> 



I v dvdVLN a N^ 



(4) 



(5) 



(6) 



where the tensor h af3 = g af3 + u a vP projects any other tensor 
into the space orthogonal to u a , namely h a ^u a = 0. In terms of 
such momen ts the radiation energy- momentum tensor T r a/3 can be 
rewritten as dHsieh & S piegel 19W 



= (E T + Vr)u a vP + F?vP + u a Ff + V T g 



(7) 



where E r and V r are the radiation en ergy density and pressure, re- 
spectively. As in iFarris et all d2008l) . we make the additional and 
strong physical assumption that the radiation is very close to be- 
ing isotropic in the comoving frame of the fluid, thus mimicking 
the conditions of the optically thick regime. However, while the ra- 
diation pressure is actually set to be V r = E T /3, as the isotropic 
assumption implies, the radiation flux is allowed to assume non- 
vanishing values, although with the constraint that F? /E T <C 1. 
Hence, the radiation field is only approximately isotropic. 

The full set of equations describing the dynamics of the sys- 
tem is 



a (3 



V Q T r 



o, 



-G 



(8) 
(9) 
(10) 



While Eqs. {8]), and J9]) represent the well known continu- 
ity equation and the energy momentum equation, Eq. dTTTb ex- 
presses the evolution of the radiation field, where G> is the ra- 
diation four-force density. The latter depends on the physical in- 
teraction between matter and radiation and is therefore specific to 
the problem considered. In full generali ty this tensor is given by 
dMihalas & Mihalasll 19841 : IShapirdl 19961) 



Gr 



dvdQ , 



(11) 



where = xt + xt an d = vt + rjt are the total opacity and 
emissivity coefficient^ each containing a thermal contribution, in- 
dicated with the superscript "t", and a scattering one, indicated with 
a superscript "s". In addition, we assume that: (i) the scattering is 
isotropic and coherent; ( ii) the thermal emissivity and the thermal 
opacity coefficients are related to the Planck function B v through 
Kirchhoff's law r\ l v — B u xl\ (Hi) that electrons and ions are main- 
tained at the same temperature; (iv) the opacity coefficients are in- 
dependent of frequency, x» — Kgp, where n g is the gray-body 
opacity. The last assumption, in particular, prevents us from taking 
into account photoionization effects, which are therefore not con- 
sidered in our analysis. 

Under these conditions, which are indeed the same considered 

3 Note that although both are referred to as "coefficients", \v and 77^ have 
different units. The dimensions of \v are cm -1 , while those of r} v are 
erg cm _3 s _1 Hz _1 sr _1 . 



by IFarris et all d2008l) « the radiation four-force can be written in 
covariant form as 



Gr = x(Er - 4nB)u a + ( x * + x)F? . 



(12) 



where 4ttB = a ra dT 4 is the equilibrium black-body intensity, with 
T the temperature of the fluid and a ra d is the radiation constant. 
The temperature is estimated from the ideal-gas EOS via the ex- 
pression 



m p p 
k B p ' 



(13) 



where, we recall, ks is the Boltzmann constant and m p the 
rest-mass of the proton . In this paper we consid er the case of 
bremsstrahlung opacity dRvbicki & L ightman 19861) 



Xbr 



1.7 x 10" 



i rri-7/2 ry2 -1 

i u Z n e rii cm 



= 1.7 X 10 -25 T -7/2^s 



(14) 



where n e and m ~ n e are respectively the number densities of 
electrons and ions (protons) expressed in cgs units, while Tr is the 
equilibrium temperature of both electrons and protons expressed in 
Kelvin. For the scattering opacity we consider the Thomson scat- 
tering opacity and we recall that the Thomson cross sections of 
electrons and protons are: ar, e = 6.6524586 x 10 _25 cm 2 and 
c^T,p = (m e /m p ) 2 <jT,e, respectively. Hence, the Thomson scatter- 
ing opacities of electrons and of protons are given by 

xt = &T,e n e = (TT,e ( — J = 0.397726 pegs cm -1 , (15) 
\m p J 

Xp = (T T ,p n p = a T , P ( — ) = 1.17968 x 10" 7 p cgs cm -1 . 

(16) 

We recall that the electron- scattering opacity domin ates over free - 
free opacity at low densities and high temperatures dHarwitlll998h . 
where the interaction between electrons and ions is weak. It is 
worth stressing that, because of the assumptions made, an incoher- 
ent process such as Compton scattering, with a cross section that is 
frequency dependent, cannot be consistently taken into account and 
it is therefore neglected. Finally, as customary, the optical thickness 
is defined as the line integral of the opacities between two points in 
the fluid 



Jo 



+ X S )ds. 



(17) 



In practice, we approximate expression (fTTt as r ~ (%* + x* )L, 
with L being a typical length scale of the problem. 

We refer to Appendix lAl for a summary about the conversion 
between cgs and geometrized units. 



2.2 Numerical methods 

We solve the equations of general relativistic non-dissipative ra- 
diation hydr odynamics ([ll-flOl thr ough a modified version of the 
ECHO code (iDel Zanna et al.ll2007l) . which adopts a 3 + 1 split of 
spacetime in which the spacetime is foliated into non-intersecting 
space-like hyper- surfaces St, defined as isosurfaces of a scalar time 
functi on t. Within this appro ach, the metric is decomposed accord- 
ing to dArnowittetal.lll962h 



ds 2 = -adt 2 + jij (dx'+/3Mt)(dx J +/3 j dt), 



(18) 
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where a is the lapse function, f5 l is the shift vector, 7^ is the spatial 
metric tensor, and 



(-a, 00, (n A 



-1), 



(19) 



is the future-pointing time-like unit vector normal to the slices E*. 
The observer mo ving with four- veloci ty n M = {1/a, — /3 1 /a} is 
called Eulerian dSmarr& Yorklll978h ." Any vector V ^ (or simi- 
larly a tensor) may be projected in its temporal component V n = 
-n^V^ and spatial component ±V M = (gt + n^n v )V v . As a 
result, any spatial vector (or tensor) must necessarily have a 
vanishing contravariant temporal component V t = 0, whereas its 
covariant temporal component is Vt = g^tV^ — PiV 1 , in general 
different from zero. The 3+1 splitting procedure just described 
can be applied to the vectors and tensor introduced so far to yield 



01 

u 




1 


(20) 


J - m 


= W a(3 + SV 


' +n a S* 3 +Un a n* 3 , 


(21) 


r 


= aF f v n a + / r a 


j 


(22) 


± r 


= Rf + Sfn p 


+ n a S I 3 + Urn !!, 13 , 


(23) 



where all the tensors v", W^ u , 5 f/x , f?, R£ v , S? correspond to 
the familiar three-dimensional quantities as measured by the Eule- 
rian observers, are purely spatial, and have indices that are raised 
and lowered by the spatial metric 7^ . In particular, the newly intro- 
duced quantities are related to the corresponding quantities in the 
comoving frame by 



D 


= pT, 


(24) 


W ij 


= phV 2 v z v* + p 7* J , 


(25) 




= phY 2 v l , 


(26) 


u 


= phV 2 - p , 


(27) 


r>ij 


= ^ r r 2 «v + r(/ r v + / r V) + v^ i3 , 


(28) 


si 


= ^ r rV + r(aF r V + /;), 


(29) 




= -E r T 2 + 2aTF r t - . 


(30) 



A few comments about the quantities in the equations above can be 
useful. The vectors v % and /* are the velocity and the radiation flux, 
respectively, as measured by the Eulerian observers, while V — 
(1 — v 2 ) -1 / 2 = au* is the Lorentz factor of the bulk flow. In 



particular, the radiation flux vector is /* = F r 2 + /3 Z F^ 
is computed from the orthogonality condition F^u a ■ 
given by 



Vifr 

a 



where F* 
= and is 



(31) 



It is interesting to note that U T — R^^n a n^ is the radiation en- 
ergy density as measured by the Eulerian observers, in analogy with 
what happens for the conserved energy density of the fluid U de- 
fined by d27l ). 

The general-relativistic radiation-hydrodynamics equations 
are then written in the following conservative form 



d t U + 3 % T l = S , 



(32) 



which is appropriate for numerical integration via standard high- 
resolution shock-capturing (HRSC) methods developed for the Eu- 
ler equations. The conservative variables and the corresponding 



fluxes in the i-direction are respectively given by 

av'D - /3 % D 



D 

U 



aW) - ftSj 
aRl j - /3i(S r )j 



, (33) 



whereas the sources, in any stationary background metric, can be 
written as 



s = 



\aW ik d jlik + SidjP* - Ud 3 a + a(G r ) 3 
\W ik p j d jlik + Wjdj? - S'dja + a 2 G\ 



Sldja 



a Cr r 



(34) 

where only purely spatial quantities are present. We note that 
y/j = y/^g/a is the determinant of the spatial metric. In our 
setup for two dimensional simulations presented in Sec. |4]we as- 
sume the metric given by the Kerr solution with the limiting case 
of Schwarzschild metric for vanishing black-hole spins. 

The radial numerical grid is discretized by choosing N r points 
from r m in to r ma x, non-uniformly distributed according to the fol- 
lowing scheme 



Ti = r min + ai tan (a 2 Xi) , 

Xi = {fi - r min )/(r max ?"min J i 



(35) 
(36) 



where a\ — (r max — r m i n )/ao, CL2 = arctanao, while ri are the 
coordinate points of the uniform grid from r m in to r m ax- In prac- 
tice, the free parameter ao controls the extent to which the grid 
points of the original uniform grid are concentrated towards r m i n , 
and we have chosen ao in the range [5 — 10] in most of our simu- 
lations. The angular grid is taken to be uniform. 

The set of hydrodynamics equations is discretized in time with 
the method of lines and the evolution is performed with a second- 
order modified Euler scheme. A fifth-order finite-difference algo- 
rithm based on an upwind monotonicity-preserving filter is em- 
ployed for spatial reconstruction of primitive variables, whereas 
a two- wave HLL Riema nn solver is used to en sure the shock- 
capturing properties (see Del Zan na et al 1 (120071) for further de- 
tails). As a final remark we note that as customary in HRSC meth- 
ods, we introduce a tenuous and static "atmosphere" in the regions 
of the fluid where the rest-mass density falls below a chosen thresh- 
ol d value. When this happens, we follow the prescription detailed 
in lBaiotti et al.l d2Q05h as far as the hydrodynamical quantities are 
concerned, while the primitive variables of the radiation field are 
frozen to the values at previous time-step. 



3 VALIDATION OF THE CODE 

The purely magnetohydrodynamic version of the code has been 
validated over the same numerical tests extensively described 
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Table 1. Description of the initial states in the shock-tube tests with radiation field. The different columns refer respectively to: the test considered, the radiation 
constant, the adiabatic index and the thermal opacity. Also reported are the rest-mass density, pressure, velocity and radiation energy density in the "left" (L) 
and "right" (R) states. 



Model 


7 


«rad 


4 


PL 


PL 




U L 


E Y ,L 


PR 


PR 


U R 


E r ,R 


1 


5/3 


1.234 x 10 10 


0.4 


1.0 


3.0 x 10- 


-5 


0.015 


1.0 x 10~ 8 


2.4 


1.61 x 10~ 4 


6.25 x lO -3 


2.51 x 10~ 7 


2 


5/3 


7.812 x 10 4 


0.2 


1.0 


4.0 x 10- 


-3 


0.25 


2.0 x 10- 5 


3.11 


0.04512 


0.0804 


3.46 x lO" 3 


3 


2 


1.543 x 1CT 7 


0.3 


1.0 


60.0 




10.0 


2.0 


8.0 


2.34 x 10 3 


1.25 


1.14 x 10 3 


4 


5/3 


1.388 x 10 8 


0.08 


1.0 


6.0 x 10" 


-3 


0.69 


0.18 


3.65 


3.59 x 10- 2 


0.189 


1.3 




Figure 1. Solution of the shock tube test 1 (left panel) and 2 (right panel) as reported in Table [T] From top to bottom the panels report the rest-mass density, 
the velocity, the radiation energy density and the radiation flux. 



in Del Zanna etal obtaining the same convergence proper- 

ties and will not be reported here for compactness. For the radiation 
part of the code, on the other hand, there are only a few analytic or 
semi-analytic tests that can be adopted, as we discuss below. 

3.1 Shock- tube problems 

Considering a flat spacetime, we have followed lFarris et all d2008h . 
who proposed and solved four shock-tube tests in which nonlin- 
ear radiation-hydrodynamic waves propagate. The initial states of 
these tests are reported in Table [T] and are chosen in such a way 
that the discontinuity front at x = remains stationary, namely it 
has zero velocity with respect to the Eulerian observer of the code. 
The values of the fluxes, not reported in Table Q] are chosen to 
be two orders of magnitude smaller than the energy density of the 
radiation field. In these tests local thermodynamic equilibrium is 
assumed at both ends x = ±X, with X = 20, and this is obtained 
by adopting a fictitious value of the radiation constant a ra d, namely 
ttrad = Et^l/Tl, which is then used to compute E y ,r = a ra d^H 
(here the indices L and R indicate the "left" and "right" states, re- 
spectively). The scattering opacity K s g is set to zero in all of the 
tests, while the value of the thermal opacity n g is reported in Ta- 
bled! 

Each test is evolved in time until stationarity is reached. The 
semi-analytic solution that is used for comparison wit h the numer- 
ical one has been obtained following the strategy by iFarris et all 



and it implies the solution of the following system of ordi- 
nary differential equations 

d x U(P) = S(P) , 

where 



( 9 \ 




/ P u* 


) 




/ 







p 




rpQx 











u x 


u = 


rpXX 




s = 













rj-\0x 








-G° 




\ Ff ) 




{ T T XX 


J 




V 


-Gf 


) 



Figures Q] and |2] show the comparison of the numerical solu- 
tion with respect to the semi-analytic one in the four cases con- 
sidered, which correspond, respectively, to the propagation of a 
nonrelativistic strong shock, of a mildly relativistic strong shock, 
of a highly relativistic wave and of a radiation pressure dominated 
mildly relativistic wave. In particular, Fig.[T|reports the solution for 
the tests 1 and 2 , which contain a true discontinuity represented by 
a shock front, while tests 3 and 4 have continuous configurations 
and are shown in Fig.|2] 

The tests have been performed with N = 800 uni- 
formly s paced grid poin t s usi ng the MP5 slope limiter de- 
scri bed inlDel Z anna et al] J2007h and a HLL Riemann solver. Un- 
like iF^ri^^n^Soi^rwe have not boosted the solution. This re- 
sults in a more stringent test for the code to maintain stationar- 
ity and it also explains why the profiles of the vector quantities, 
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namely t he velocity and th e radiation flux, do not match those 
shown bv lFarris et al.l d2008l) . The numerical solution is almost in- 
distinguishable from the semi-analytic one in all of the profiles re- 
ported in the figures, thus proving the ability of the code in handling 
different physical regimes of the radiation field within an optically- 
thick approximation. 



4 BONDI-HOYLE ACCRETION FLOWS 
4.1 Initial and boundary conditions 

Our attention is focused on a Bondi-Hoyle accretion flow onto a 
black hole of galactic size with M B h = 3.6 x 10 6 M©, that we 
investigate by performing numerical simulations on the equatorial 
plane, i.e. = tt/2. Despite the l ong history i n literature on this 
type of accretion (see the review by Edgar (2004)), no stationary so- 
lution for a radiation-hydrodynamics Bondi-Hoyle flow is known, 
and which could have been used as suitable initial data. As a result, 
we let the code converge to the nearest stationary solution after 
specifying the hydrodynamical solution of the Bondi-Hoyle flow, 
to which we add a radiation field with uniform and small energy 
density E Y . 

Most of our discussion hereafter refers to accretion onto 
Schwarzschild black holes, although also rotating black holes will 
be briefly presented in Sect. 14.3.31 The code solves therefore the 
equations in a general Kerr metric expressed in Boyer-Lindquist 
coordinates, so that the initial velocity field, specified in ter ms of 
an asymptotic velocity Voo, is given by jFont & Ibanez 1998) 

v r — y/Y^Voo COS (j) , (37) 

v* = -v / 7^^oosin0. (38) 

These relations guarantee that the velocity of the injected gas 
at infinity is parallel to the x— direction, while v 2 = ViV 1 = v 2 ^ 
everywhere in the flow. Other quantities that need to be set initially 
are: the asymptotic sound speed c S)0 o, and the asymptotic pressure, 
from which the asymptotic rest-mass density poo follows directly 
(see values reported in Table|2]for all of the models considered). For 



all the simulations we will consider a gas of nonrelativistic elec- 
trons and hence with an adiabatic index 7 = 5/3. The velocities 
used in our models and presented in Table |2] are chosen to be suf- 
ficiently high so as to open a shock cone (see details below). Any 
chosen Voo implies a restricted range of asymptotic sound speeds, 
if a reasonable Mach number should be considered. We remark that 
our models do not aim at modelling any specific astrophysical sce- 
nario, but rather at highlighting the role of the back-reaction of the 
radiation in an optically thick, relativistic Bondi-Hoyle accretion 
flow. 

Similarly, the radiation field is initialized to a value such that 
the radiation temperature T rad = (E Y /a ra d) 1/4 ~ 1.5 x 10 5 K. 
While this may seem an arbitrary choice, we have verified through 
a series of numerical simulations that, on long-term evolutions, the 
value of the obtained luminosity is not dependent of this initial 
choice. The computational grid consists of N r x N<f, numerical cells 
in the radial and angular directions, respectively, covering a compu- 
tational domain extending from r m i n = 2.1 M to r max = 200 M 
and from m i n = to ma x = 2tt. For our fiducial simulation we 
have chosen N r = 1536 and N<f, — 300, but have also verified that 
the results are not sensitive to the resolution used or to the location 
of the outer boundary. 

The boundary conditions in the radial direction are such that 
at the inner radial grid point we implement inflow boundary con- 
ditions by a simple zeroth-order extrapolation {i.e. a direct copy) 
of all variables. At the outer radial boundary, on the other hand, we 
must distinguish between the upstream region {i.e. with tt/2 < (j) < 
3/2?r), and the downstream region {i.e. with — ty/2 < (j) < tt/2). 
In the upstream region we continuously inject matter with the ini- 
tial velocity field of d37])-([38t, thus reproducing a continuous wind 
at large distances, while in the downstream region we use outflow 
boundary conditions. Finally, symmetric {i.e. periodic) boundary 
conditions are adopted at = 0. The simulations are performed 
with a Courant-Friedrichs-Lewy coefficient that may vary accord- 
ing to the model and it is typically in the range ~ [0.01, 0.5]. 

In addition to the "classical" Bondi-Hoyle initial data, we will 
also consider a set of simulations in which the thermodynamics of 
the flow is slightly altered in order to reduce the temperature of the 



GR radiation hydrodynamics of accretion flows 7 



Table 2. Initial models adopted in numerical simulation. From left to right the columns report: the name of the model, the asymptotic flow velocity Vqo, 
the asymptotic sound speed c s ,oo, the asymptotic Mach number .Moo, the initial temperature, the initial rest-mass density and the accretion radius r a = 
GM / (v^ + OQ ). Perturbed Bondi-Hoyle solutions are generated through injection of low pressure gas as described in the text. Each model is evolved until 
stationarity is reached, and in any case up until at least t = 20000 M. The adiabatic index was set to 7 = 5/3 and the black hole spin a = 0. The mass of the 
black hole is M BH = 3.6 x 10 6 M© and the radial grid extends from r min = 2.1 M to r max = 200 M. 



Model 


^00 


Cs, 00 


Moo 


T[K] 


Poo 1 


cgs 


] 


r a [M] 


V08.CS07 


0.08 


0.07 


1.14 


3.22 


X 


10 10 


3.22 


X 


10" 


-12 


88.5 


V09.CS07 


0.09 


0.07 


1.28 


3.22 


X 


10 10 


3.22 


X 


10" 


-12 


76.9 


V10.CS07 


0.10 


0.07 


1.42 


3.22 


X 


10 10 


3.22 


X 


10" 


-12 


67.1 


V11.CS07 


0.11 


0.07 


1.57 


3.22 


X 


10 10 


3.22 


X 


10" 


-12 


58.8 


V07.CS06 


0.07 


0.06 


1.16 


2.36 


X 


10 10 


4.39 


X 


10" 


-12 


117.6 


V07.CS07 


0.07 


0.07 


1.0 


3.22 


X 


10 10 


3.22 


X 


10" 


-12 


102.0 


V07.CS08 


0.07 


0.08 


0.77 


4.22 


X 


10 10 


2.45 


X 


10" 


-12 


88.5 


V07.CS09 


0.07 


0.09 


0.77 


5.35 


X 


10 10 


1.93 


X 


10" 


-12 


76.9 


p.V09.CS07 


0.09 


0.07 


1.28 


3.22 


X 


10 9 


3.22 


X 


10" 


-12 


76.9 


p.V10.CS07 


0.10 


0.07 


1.42 


3.22 


X 


10 9 


3.22 


X 


10" 


-12 


67.1 


p.Vll.CS07 


0.11 


0.07 


1.57 


3.22 


X 


10 9 


3.22 


X 


10" 


-12 


58.8 


p.V18.CS07 


0.18 


0.07 


2.57 


3.22 


X 


10 9 


3.22 


X 


10" 


-12 


26.8 



gas. We denote these models as "p-models" in Table[2] In essence, 
the perturbed Bondi-Hoyle accretion flows are obtained by inject- 
ing gas of lower pressure than required by the stationary solution 
at the upwind boundary, with a proportionally reduced radiation 
energy density E T (see Sect. 14.3.21 for details). 



4.2 Computation of luminosity 

The key new quantity that our code allows to compute is the emit- 
ted luminosity. Since the code explicitly calculates the radiation 
fluxes fl at each time step, we use them to compute the intrinsic 
luminosity emitted from the optically thick region as 

L= / flriid S pt , (39) 
Jn 

where S op t is the surface of the volume Q enclosing an optically 
thick region within the computational domain, while the scalar 
product fl ra provides the projection of the local radiation flux onto 
the normal to the radiating surface. Because of the nearly isotropic 
assumption made for the radiation field and because of the rough 
spherical symmetry of the physical system under consideration, the 
fluxes in the angular directions are expected to to be much smaller 
than the radial ones and to almost cancel. As a result, and for sim- 
plicity, we approximate the scalar product above as /* ra — / r r , 
thus computing the luminosity as 

L = 2 E Ur) n A^n] \r=l , (40) 

n=l 

where A0 n is the angular size of a grid cell and we perform the sur- 
face integral at the radial position of the last optically-thick surface, 
i.e. where r = 1; the factor 2 accounts for both the contributions 
above and below the equatorial plane. 

The luminosity computed in this way comprises two different 
contributions. The first one is an accretion-powered luminosity that 
it is directly proportional to the mass-accretion rate M through a 
relation of the type L acc = rjM, where the coefficient 77 expresses 
the efficiency of the conversion of gravitational binding energy into 
radiation. The main dissipative mechanism is provided by compres- 



sion of the fluid when this has nonzero thermal conductivity A 
second contribution to the total luminosity (|40l is given by dissipa- 
tive processes related to shock heating that, as we will show below, 
can provide a considerable contribution to the total emission. 

However, since we are dealing with inviscid non-magnetized 
fluids, the luminosity (Rot obviously cannot provide the contri- 
bution coming from dissipative processes driven by viscosity (of 
whatever origin), and that can be a significant part of the accretion- 
powered luminosity in a realistic accretion scenario. We recall, for 
instance, that in the classical Shakura-Sunyaev thin-disc model the 
main dissipative mechanism comes from the viscous stress ten- 
sor, directly proportional to the total pressure via the "alpha" pa- 
rameter. Similarly, in spherical accretion, a realistic viscous fluid 
with nonzero bulk viscosity will produce a viscous dissipation 
adding to the one coming from the fluid compression. In sum- 
mary: in realistic accretion scenarios one should expect that both 
thermal conductivity and viscosity act as transport coefficients of 
dissipative processes and lead to contributions to the emitted lu- 
minosity. In our treatment, however, only the effects of the for- 
mer one can be accounted for. Hereafter, the luminosities and the 
accretion rates will be reported in Eddington units, i.e. LEdd = 
4irGMm p c/(TT,e ^ 1.26 x 10 38 (M/M ) erg s~\ M E dd 
L E dd/c 2 ~ 1.39 x 10 17 (M/Mq) gs" 1 . See also ggj for the 
Eddington luminosity in the geometrized units of the code. 



4.3 Results 

Before entering into the details of our results, it is useful to briefly 
review the main features of the relativistic Bondi-Hoyle accre- 
tio n as investigated th r ough purely hydrodyn a mical simulations 
by IPetrich et al l dl989h : iFont & Ibanezl dl998h : iFont et all dl998h 
and lFont et al.l (1 19991) . Overall, these studies have highlighted that 
when a homogeneous flow of matter moves non-radially towards a 



4 We recall that the thermal conductivity is related to the opacity and its 
effects are therefore accounted for in our analysis. For instance, the thermal 
conductivity computed using the ordinary diffusion appro ximatio n of stellar 
interiors is given by xt = (4/3)a rad cT 3 /x s ( Schwartz ! 1967b . 
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Figure 3. Rest-mass density in cgs units on a logarithmic scale for model V09.CS07 in a purely hydrodynamical evolution (left panels) and in a radiation- 
hydrodynamics evolution (right panels). Different rows refer to different times of the evolution and white regions correspond to densities slightly below the 
threshold for the colour code at around 10 -12 g/cm -3 . Note that the presence of a radiation field reduces the rest- mass density considerably near the black 
hole, suppressing the accretion rate. 
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Figure 4. Le/j? Panel: logarithm of the ratio of radiation pressure over gas pressure for the model V09.CS07 at early times. Right Panel: the same as the right 
panel but at later times, when stationarity had been reached. 



compact object, a shock wave will form close to the accretor. De- 
pending on the adiabatic index and on the asymptotic Mach num- 
Moo, the shock can either come very close to the accretor 
or be a t a certain distance from it [see, for instance, iFoglizzo et al.1 
J2005h l. In general, for any given value of the adiabatic index, there 
is a minimum asymptotic Mach number above which a shock wave 
of conic shape, i.e. a "shock cone", forms downstream of the accre- 
tor. On the other hand, asymptotic Mach numbers below the criti- 
cal value produce a shock wave that, initially formed in the down- 
stream region, opens progressively and reverses in the upstream re- 
gion as a bow shock. More recently, two different studies have shed 
additional light on the p hysics of relati v istic B ondi-Hoyle accretion 
flows. In the first one, iDonmez et al.l d201Qh . reported the occur- 
rence of the so called flip-flop instability of the shock cone in the 
relativistic regime and have also shown that quasi-periodic oscilla- 
tions ofsonic nature are produced in the shock cone. In the second 
one, IPennerl ( 1201 lh investigated the effects of a uniform magnetic 
field, finding that it produces an increase in the cone opening angle 
and in the mass accretion rate. 



4. 3. 1 Classical Bondi-Hoyle accretion 

We start our analysis by considering the extent to which a radiation 
field affects the dynamics of the classical Bondi-Hoyle flow, com- 
paring the dynamics for very similar physical conditions. The ini- 
tial models, which are the first seven reported in Table [2] have very 
high temperatures and, consequently, high thermal conductivitieo 
As mentioned above, for any given value of the adiabatic index, 
there is a critical asymptotic Mach number «A4oo, c , usually close to 



5 We recall that the relativistic Mach number is defined as 
M = Tv/(c s T s ), where T and T s are the Lorentz factors of the 
flow and of the sound speed, respectively. 

6 The present version of the code does not allow to handle stiff source terms 
that aris e in the r adiation-hydrodynamics equations when the conductivity 
is small (Szoke et al. 2006). Work is in progress to cope with this difficulty. 



unity, above which a shock cone forms in the downstream region 
and below which the shock cone reverses in the upstream region. 
Our simulations indicate that, for values of the Mach number close 
to the critical one, the radiation effects on the dynamics are most 
evident. This is shown in Fig. [3] for model V09.CS07, where we 
have reported the distribution of the rest-mass density at three dif- 
ferent times in a purely hydrodynamical evolution (left panels) and 
in a radiation-hydrodynamic evolution (right panels). This model, 
in particular, provides an example in which the radiation field pre- 
vents the reversal of the shock cone from the downstream region 
into the upstream region, which instead takes place in the purely hy- 
drodynamical evolution. Since the dynamics of V09.CS07 becomes 
radiation-pressure dominated around t ~ 5000 M, the explanation 
of this effect is simple: In such conditions the effective adiabatic 
index of the fluid-plus-radiati on medium is smaller than that of the 
fluid alone [see Eq. (70.22) of iMihalas & Mihalasl dl984h l 



7eff 



5/2 + 20g + 16g 2 
(3/2 + 12g) (! + <?) ' 



(41) 



where q = V T /p. This fact has two important consequences. The 
first one, which we will discuss shortly when commenting Fig. [7] 
is to increase the rest-mass density jumps across shock fronts. The 
second one, is exactly to favour the generation of the shock cone 
downstream of the accretor, as firstly not iced bv lRuffertldl996h and 
later confirmed bv lFont & Ibanezld 19981) . 

As clearly shown in Fig. [3] the radiation-hydrodynamics evo- 
lution of model V09.CS07 is remarkably different from the purely 
hydrodynamical one, and it can be divided in the following stages. 
After the shock cone has fully opened in the downstream region 
(top right panel of Fig. [3]), the flow becomes radiation-pressure 
dominated, making the shock cone oscillate from one side of the 
accretor to the other, in a way that resembles the flip-flop in- 
stability already enco untered in relativistic Bondi-Hoyle flows by 
IDonmez et all bold) . This transient behaviour, captured in the 
right-middle panel of Fig. [3] is accompanied by an outflow of mat- 
ter expelled by radiation pressure beyond the computational grid. 
After that, the system relaxes to a stationary configuration charac- 
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Figure 5. Rest-mass density in cgs units on a logarithmic scale for model V10.CS07 in a purely hydrodynamical evolution (left panels) and in a radiation- 
hydrodynamics evolution (right panels). Different rows refer to different times of the evolution and white regions correspond to densities slightly below the 
threshold for the colour code at around 10 -12 g/cm -3 . Note that the presence of a radiation field reduces the rest- mass density considerably near the black 
hole, suppressing the accretion rate. 
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Figure 6. Evolution of the mass accretion rate in Eddington units for models V07.CS07 and V07.CS09 (left panels) and V09.CS07 and V10.CS07 (right panels). 



terized by the presence of a shock cone with a much smaller open- 
ing angle than in the hydrodynamical solution, giving rise to a "re- 
duced" shock cone. Also shown in Fig. [4] is the ratio between the 
radiation pressure and the fluid pressure, reported in the left panel 
for an early and fluid-pressure dominated stage of the evolution, 
and in the right panel for a late and radiation-pressure dominated 
one. 

A very similar behaviour to the one discussed so far is shown 
in Fig.|5]for model V10.CS07, the initial Mach number of which is 
only slightly larger than model V09.CS07. However, in this case the 
higher fluid velocity causes supercritical behaviour both in the hy- 
drodynamical and the radiation-hydrodynamical evolution so that 
the shock cone remains in the downstream region. The close sim- 
ilarity between the dynamics of models V09.CS07 and V10.CS07 
in the presence of the radiation field is also testified by the asymp- 
totic mass accretion rate, which is M ~ 13.14MEdd for model 
V09.CS07 and M ~ 10.24M E dd for V10.CS07. 

An information complementary to that of Fig. [3] - [5] is pro- 
vided by Fig. [6] which shows the evolution of the mass accre- 
tion rate for a few selected models. For each of these models 
both the purely hydrodynamical evolution (red solid lines) and the 
radiation-hydrodynamical one (blue dashed lines) are considered. 
A few comments are worth making about this figure. The first one 
is that, once stationarity is reached, the mass accretion rates of 
the radiation-hydrodynamics models are significantly smaller than 
those of the corresponding hydrodynamics models. This result was 
of course expected, because of the obstructive effect of the radiation 
pressure. The second comment is that the reversal of the shock cone 
in the hydrodynamics models V09.CS07, V07.CS07, V07.CS09 and 
in the radiation-hydrodynamics models V07.CS07 and V07.CS09 
leads to an increase of M, as highlighted by the arrows. For the 
hydrodynamical version of model V09.CS07, for instance, this in- 
crease starts at t ~ 12000 M, as reported in the top-left panel of 
Fig. [6] Finally, we find that all models accrete at super-Eddington 
rates even when a radiation field is present. This is not surprising, 
since the Eddington limit holds strictly only in spherical symmetry, 
which is not fulfilled in wind- like accretion. Moreover, it should be 
remarked that the classical Eddington limit is computed in a frame- 




Figure 7. Comparison of the rest-mass density jump across the shock front 
at time t = 9500 for models V09.CS07 in a purely hydrodynamics evolu- 
tion (red solid line) and in a radiation hydrodynamics evolution (blue dashed 
line). The location of the shock was re-normalized to lie at x = and is 
displaced by Ax ~ 0.2M between the two runs. 



work where only the electron Thomson cross section contributes to 
the radiation pressure. 

Additional differences between the hydrodynamics and the 
radiation-hydrodynamics evolutions emerge after comparing the 
jumps experienced by the rest-mass density across a shock wave 
in a representative model. Such a comparison is reported in Fig. [7] 
for model V09.CS07, showing the variation of the rest-mass density 
across the shock that is produced at time t — 5000M and visible 
in the two top panels of Fig. [3] The two curves have been obtained 
after slicing the rest-mass density along an "x-direction" perpen- 
dicular to the shock front, and sliding the two profiles so that the 
shock is located at the same x = position for both the hydrody- 
namical and the radiation-hydrodynamical evolution. Negative and 
positive values of the x-coordinate refer therefore to the unshocked 
and to the shocked region, respectively, while the rest-mass density 
has been normalized to the value in the unshocked region. 

The first comment about this figure is that the density jump in 
the hydrodynamics evolution is slightly smaller than the value of 4 
predicted by the theoretical expectation of P2/P1 ~ (7+ — 1) 
valid for an ideal gas EOS. This effect may be due to the presence of 
both numerical diffusion and tangential velocities along the shock 
front. The second comment is that the compression ratio across the 
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Figure 8. Left Panel: logarithm of the optical thickness for the model V09.CS07 once stationarity is reached. Right Panel: Logarithm of the modulus of the 
radiative flux (in geometrized units) in the same model at the same time. 



Table 3. Mass accretion rate, luminosity and efficiency r\ Bri as defined in 
Eq. J42l of Bondi-Hoyle accretion in the quasi- stationary regime. 



Model 


Moo 


Macc/M E dd 


£/^Edd 




V08.CS07 


1.14 


17.9 


3.10 


0.14 


V09.CS07 


1.28 


13.14 


4.11 


0.22 


V10.CS07 


1.42 


10.24 


6.18 


0.36 


V11.CS07 


1.57 


8.32 


6.77 


0.38 


V07.CS07 


1.0 


27.64 


1.44 


0.05 


V07.CS08 


0.87 


21.74 


2.14 


0.09 


V07.CS09 


0.77 


17.44 


1.91 


0.10 


p.V10.CS07 


1.42 


11.58 


5.35 


0.29 


p.Vll.CS07 


1.57 


8.34 


7.05 


0.40 


p.V18.CS07 


2.57 


4.02 


30.5 


0.69 



shock increases by 3% in the transition from a hydrodynamical 
to a radiation-hydrodynamical evolution. This effect can again be 
understood by regarding the fluid in the radiation-hydrodynamics 
evolution as an effective fluid having a smaller adiabatic index, as 
indeed expected in the radiation-pressure dominated regim e. This 
result is al so in agreement with analy tical investigations by iGuessI 
dl96Qh and lMihalas & MihalasUl984) ($104). 

Before computing the luminosity as described in Sec. 14.21 it is 
important to make sure that the physical conditions chosen corre- 
spond to those required by the code, namely the presence (and the 
persistence) of an optically thick regime. Of course all of the mod- 
els considered in our simulations and reported in Tableware in such 
a physical regime, with only very limited regions where the optical 
thickness can be ~ 0(1) during the evolution. As a representative 
example, the left panel of Fig. [8] shows the optical thickness when 
the system as relaxed to stationarity, for the same model V09.CS07 
that we have extensively described so far. The right panel of Fig. [8] 
on the other hand, shows the corresponding intensity of the momen- 
tum of the radiation field. After comparing with the right-bottom 



panel of Fig.[3j it is easy to realize that the distribution of the radia- 
tive fluxes is obviously correlated with the rest-mass density dis- 
tribution, but also that a good portion of the radiative emission is 
concentrated along the shock fronts of the reduced shock cone. 

The evolution of the emitted luminosity and of the mass- 
accretion rates are illustrated in the two panels of Fig. [9] More 
specifically, the left panel, which reports models with increasing 
Mach number but having the same initial temperature, shows that 
the luminosity increases with Moo and reaches stationary values 
of a few Eddington units. On the other hand, the right panel, which 
reports models with the same asymptotic velocity but different tem- 
peratures, shows that stationarity is reached on longer timescales 
and a correlation with the final luminosity is less robust. In all of 
the models shown, the first bump around t ~ 2000 M is due to the 
initial opening of the shock cone. 

By providing the first self-consistent computation of the lumi- 
nosity in a Bondi-Hoyle accretion flow, our calculations allow us to 
derive the efficiency of the accretion flow r] Bn . We remark that the 
concept of rj BrL for a Bondi-Hoyle flow, with a nonzero velocity of 
the matter at infinity, is not the same as in standard accretion discs, 
where the gas flow is supposed to start from matter at rest at infin- 
ity. Thus, we define an effective Bondi-Hoyle luminosity efficiency 
Vbu as 

ri^ = — , (42) 

M acc C 2 + ±Moo«§o 

where the denominator takes into account a kinetic contribution to 
the energy flux. We report the values of M &cc / M^dd, L/L^dd and 
rj BU in Table [3] for those models presenting a quasi- stationary ac- 
cretion pattern. From the data reported in the Table it is possible to 
deduce the existence of two different regimes in a radiative Bondi- 
Hoyle accretion flow. A first regime, corresponding to Moo < 1, 
where the luminosity is dominated by the accretion-powered lumi- 
nosity and thus proportional to M. A second regime, correspond- 
ing to Moo > 1, where the luminosity is instead dominated by 
the emission at the shock front. In particular, by comparing the first 
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Figure 9. Luminosity and mass-accretion rates in Eddington units in classical Bondi-Hoyle accretion flows. The left panel collects models with different initial 
velocities but with the same sound speed. In contrast, the right panel collects models with the same initial velocities but with different values for the sound 
speed. All simulations were evolved until they reached stationarity, and up to t = 30000 M at most. 



four models that have the same initial asymptotic sound speed, we 
note that, as the asymptotic Mach number is increased, the accre- 
tion rates decrease. This effect is due to the reduced opening angle 
of the shock cone. The corresponding luminosity, on the other hand, 
increases, because of the enhanced dissipation at the shock front. 

As a final remark we note that, as already discussed in Sec. 14.21 
the luminosities we have reported here can only provide lower lim- 
its on the energy efficiency rj Bn . We have in fact neglected not only 
viscous dissipative processes from the accretion flow, but also any 
nonthermal emission, such as inverse Compton or synchrotron ra- 
diation, which could arise from a corona developing near the black 
hole. 



4. 3. 2 Perturbed Bondi-Hoyle flow 

As mentioned in Sect. 14.11 in addition to the standard and station- 
ary Bondi-Hoyle flows, we have also considered initial conditions 
that would lead to perturbed Bondi-Hoyle accretion patterns (we 
recall that we have tagged these as the "p-models"). The rationale 
behind this choice is that of investigating how the accretion flows 
varies when the initial conditions are no longer those ensuring a 
stationary flow. At the same time, this allows us to consider mod- 
els that have lower temperatures and, consequently, lower thermal 
conductivities. 

In practice, we trigger the perturbation of the Bondi-Hoyle 
flow by acting on the thermodynamic conditions of the fluid in the 
upstream region and by producing models with values of the initial 
temperature that are typically one order of magnitude smaller than 
those in standard Bondi-Hoyle models. Because of the perturba- 
tion introduced, the dynamics of the perturbed models is typically 
characterized by a very dynamical phase before quasi- stationarity is 
reached. However, in spite of these violent transients, the perturba- 
tion introduced does not destroy the general Bondi-Hoyle pattern, 
which is recovered eventually. 

Among the perturbed models, p.V09.CS07 has the minimum 
Mach number, and is also the only one producing a shock cone that 
progressively reverses into the upstream region as a bow shock. 



The remaining three models, which all have higher Mach numbers, 
develop the usual shock cone downstream of the black hole. This 
behaviour is reported in the four panels of Fig. [10] showing the 
evolution at different times of the rest-mass density for the model 
p.V10.CS07 in a radiation-hydrodynamical evolution. Note that the 
accretion cone, that is fully formed at time t ~ 3000 M, is highly 
unstable and it goes through a rapid sequence of oscillations gener- 
ating an undulated stream in the wake. Finally, the system reaches a 
quasi-equilibrium state characterized by a reduced shock cone sim- 
ilar to that already encountered in the dynamics of standard models. 

The extraction of the light-curve and the computation of all re- 
maining quantities follows the same procedure used in the standard 
models and we have reported the mass-accretion rates and the light- 
curves in the two panels of Fig. QT] Note that the general features 
in the light-curves for the standard models are also present for the 
perturbed models. In particular, there is an initial rise in luminosity 
which corresponds to the formation of the shock cone. After that, 
between t ~ 1000 M and t ~ 2000 M depending on the model, 
a peak is produced in the light-curve which is due to the shock 
cone changing its geometry to an open cone. Interestingly, even the 
efficiency rj Bn of the perturbed models are very similar to those 
of the corresponding standard ones. For the model p.Vll.CS07, 
for instance, rj Bn = 0.40, to be compared with rj Bn = 0.38 of 
V11.CS07. 

We remark that the fluid temperature within 25 M from the 
black hole decreases more rapidly for high Mach numbers, so that 
the build up of the radiation pressure is faster for the highest Mach 
number. It should also be noted that while all perturbed mod- 
els are radiation-pressure dominated in the upstream region after 
t rsj 10000 M, this regime is reached at different times by different 
models. Furthermore, even when radiation pressure dominates the 
dynamics, there could be isolated portions of the flow where the gas 
pressure is not completely negligible. This is the case, for instance, 
in the undulated downstream part of the flow, where the ratio of gas 
pressure to radiation pressures ratio can be as high as p/V Y ~ 0.1. 

The dominant role played by the radiation pressure is im- 
printed on the accretion rate for the p-model p.V18.CS07, as 
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Figure 10. Rest-mass density in cgs units on a logarithmic scale for the perturbed Bondi-Hoyle model p.V10.CS07 at four different times in a radiation- 
hydrodynamics evolution. The two rows refer to different times of the evolution and white regions correspond to densities slightly below the threshold for the 
colour code at around 10 -12 g/cm -3 . Note that a highly dynamical transient precedes the development of a stationary flow. 



it is clear from Fig. QT| This model features the lowest quasi- 
equilibrium accretion rate and the highest luminosity. In general, 
we have found that the higher the Mach number, the higher the 
radiation pressure, and the smaller the average density around the 
black hole. The perturbed model p.V18.CS07 shown in Fig.[lOj for 
instance, has a rest-mass density which is a factor 100 smaller than 
that in model p.V09.CS07, which has the minimum Mach num- 
ber among the perturbed models and the longest relaxation time 
(cf. Fig. QT|). At the same time, the accretion rate of p.V09.CS07 
does not show the typical decline up until t = 20000 M, although 
it is radiation-pressure dominated everywhere in the numerical do- 
main. It is possible that the behaviour of model p.V09.CS07 would 
change, with the mass-accretion rate decreasing and the luminosity 
increasing, if the evolution was carried on a much longer timescale. 



4.3.3 Spinning black holes 

Although the results presented so far refer to Schwarzschild black 
holes, a number of different simulations have been performed also 
for spinning black holes, with dimensionless spin parameters rang- 
ing between and 0.999. The interest, in these cases, was that of 
determining the influence that the black-hole spin may have on the 
flow pattern and on the emission properties, for both the classical 
Bondi-Hoyle configurations and the perturbed ones. 

Overall, the modifications introduced by the black-hole spin 
are not particularly large to deserve a dedicated discussion. More 
specifically, as far as the dynamics is concerned, we have con- 
firmed that as the spin of the black hole is increased, the shock 
cone that may form in the downstream part of the fl ow is pro- 
gressively wrapped (this was originally pointed out by dFont et al. 
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Figure 11. Left Panel: Luminosity in Eddington units for the perturbed Bondi-Hoyle accretion flows. Right Panel: Mass accretion rates in Eddington units for 
the same models in the left panel. 



This distortion, however, is evident only in the immediate 
vicinity of the horizon, and typically below r ^ 20 M. Further- 
more, no significant change has been found, either qualitatively or 
quantitatively, in the luminosity, to the point that the light-curves 
for different black hole spins overlap to within 1%. These results 
suggest that if spin-related signatures in the electromagnetic emis- 
sion should exist and can be extracted, these will become evident 
only when a more sophisticated modelling of the emission pro- 
cesses {e.g. through inverse Compton in a rarefied corona) will be 
considered. This will be part of our future work. 



4.3.4 Impact on electromagnetic counterparts of supermassive 
black-hole binaries 

Considerable attention has been recently dedicated to the possi- 
bility of detecting the electromagnetic counterpart of inspiral and 
merger of supermassive binary black holes (SMBBHs) systems. 
Such a detection would not only confirm the gravitational- wave de- 
tection and help localize the source on the sky, but it would also 
provide a new tool for addressing a number of astrophysical ques- 
tions (see, e.g. lHaiman et al ]|200i)). These include the possibility 
of testing models of galaxy mergers and cl ues on the mass dist ribu- 
tion of supermassive black holes (see, e.g. ISesana et~al ] (1201 ll) and 
references therein). 

Computing the EM counterpart to the inspiral of such a bi- 
nary is not an aspect of our investigation and the physical con- 
ditions considered here badly match those expected in realistic 
scenarios describing this process, to which we plan to dedicate 
a separate investigation. The results obtained here, however, can 
serve to shed some light about a common approximation made 
in numerical simulations aim ed at estimating the luminosity from 
binary black hole mergers dBode et"aD 2010t iFarris et all 20 1C 



20091 : IZanotti et alJ 



2010 



O'Neill etal. 2009; Megevand et al 

Bode et aDl201ll ; lFarris et al.ll201ll) . All these works computed the 
bremsstrahlung luminosity without taking the back-reaction of the 
radiation into account, but rather performing a volume integral of 
the bremsstrahlung emissivity. First efforts to improve this treat- 
ment, but still without a proper radiation transfer, were initiated 



in Newtonian physics by Corral es et al ] d2010h by enforcing an 
isothermal evolution. 

To prove our conjecture that the estimates made so far in terms 
of the bremsstrahlung luminosity are optimistic, providing cooling 
times that are too short, we have computed the bremsstrahlung- 
luminosity emitted in the classical Bondi-Hoyle accretion of model 
V09.CS07 following the general relativistic prescription adopted by 
the works cited above, namely 



LBR^xlO^K^VrVTdv)^) e rg/ s > (43) 



and compared the results obtained using the estimate (143] with 
those obtained through our radiative-transfer treatment. In addi- 
tion, we have also considered an alternative calculation in which 
an isothermal evolution is enforced, and where it is assumed that 
all the changes in the temperature that are due to a local com- 
pression are dissipat ed as radiation This idea, proposed in New- 
tonian framework bv lCorrales et al.l (120101), has bee n extended to a 
general-relativistic context by IZanotti et al.l d2010L and used also 
here for comparison. 

The result of this comparison is shown in Fig|T2] where we 
have reported the three light-curves computed according to the 
approaches just described. When stationarity is reached, i.e. at 
t = 20000 M we find that L B R/£ Ed d = 78. This number 
should be contrasted with the result obtained through our self- 
consistent radiation-hydrodynamics simulations, which instead in- 
dicate L/LEdd = 4.11 (cf. Table [3])- Interestingly, the luminosity 
obtained through the isothermal approximation provides a much 
smaller value, i.e. L/L^dd = 0.09. 

While this analysis is not exhaustive and has been performed 
in the specific scenario of an optically thick Bondi-Hoyle accretion, 
it does point out that the predictions made using the simplistic es- 
timate of the bremsstrahlung luminosity via Eq. (|43] provide light- 
curves that are a factor ~ 20 larger than those obtained with a more 
rigorous approach. 
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Figure 12. Comparison among light-curves computed with different 
approaches. Solid red: luminosity obtained with the full radiation- 
hydrodynamics evolution according to Eq. ( 140) . Dashed blue: luminosity 
obtained from Eq. d43l Long-dashed green: luminosity obtained through 
the isothermal evolution approximation. See text for more explanations. 



5 CONCLUSIONS 

We h ave implemented and solved in an extension of the ECHO 
code toel Zanna et al]|2007h the equations of relativistic radiation 
hydrodynamics in the optically thick regime and on a fixed black- 
hole s pacetime when th ese equations are written in a conservation 
form (Farris et al .120081) . Within a 3 + 1 split of spacetime, we have 
discretized in time the set of equations with the method of lines 
and performed the evolution in time with a second-order modified 
Euler scheme. A fifth-order finite-difference algorithm based on 
an upwind monotonicity-preserving filter was employed for spa- 
tial reconstruction of primitive variables, whereas a two- wave HLL 
Riemann solver was used to ensure the shock-capturing properties. 
The new scheme has been successfully validated through a series 
of tests involving radiative shock tubes. 

As a first application of the new code we have considered 
the emission properties of a hot Bondi-Hoyle accretion flow onto 
a black hole with the opacity given by Thomson scattering and 
thermal bremsstrahlung only. By considering different models with 
initial temperatures around T ~ 10 10 if, an ideal-gas EOS with 
adiabatic index 7 = 5/3, and various sub-sonic and super-sonic 
regimes, we have found that the inclusion of radiation drastically 
alters the well known dynamics of Bondi-Hoyle flows in all mod- 
els considered. In particular, the system quickly enters a radiation- 
pressure dominated regime, characterized by mass accretion rates 
that, once stationarity is reached, decrease by one or two orders 
of magnitude with respect to the purely hydrodynamical evolution. 
Nevertheless, the measured accretion rates are found to be always 
super-Eddington and as high as M/M^dd ~ 25. This is in agree- 
ment with the expectation that the Eddington limit should hold 
strictly only in spherically symmetric flows. In addition, because 
the effective adiabatic index in the radiation dominated pressure 
regime is smaller than the nominal one of the gas, the radiation can 
prevent the reversal of the shock cone that is typical of Bondi-Hoyle 
flows with low Mach numbers. 

By computing the emitted luminosity through a surface in- 



tegral over the radiative fluxes at the last optically thick surface, 
our approach has allowed the first self-consistent computation of 
the light-curves for the Bondi-Hoyle flow, finding luminosities 
L/L^dd — 1 — 7. These results have been found to be independent 
of the initial conditions chosen for the intensity of the radiation en- 
ergy density. 

In addition to the classical Bondi-Hoyle accretion flows, we 
also performed simulations with perturbed setups, by injecting 
lower-temperature matter in the upstream region of the flow, which 
lead to hi ghly-dynamical trans ients reminiscent of the flip-flop 
instability (Foglizzoet "all l2005h . Although the qualitative evolu- 
tion of the accretion flow remains unchanged, the decreased initial 
temperature increases the timescale over which the flow becomes 
radiation-pressure dominated and the accretion settles in a quasi- 
stationary state. In spite of these differences, we have found that 
the main features of the Bondi-Hoyle solution, such as the pres- 
ence of the shock cone, persist under a wider class of physical con- 
ditions, even in situations departing from stationarity. Overall, our 
re sults confirm and e xtend related Newtonian studies, such as those 
bv lKlevetal.Ul995h . 

Since we have shown that the luminosity is critically affected 
by the evolution of the coupled system of hydrodynamic and ra- 
diation equations, significant changes in the luminosities should 
be expected in those scenarios which have so far been modeled 
through a-posteriori calculations to a purely hydrodynamical evolu- 
tion. A first example in this respect is given by multi-colour black- 
body spectra, while a second example is represented by the cal- 
culation of electromagnetic counte rpart to the inspiral of super- 
massive binary black-hole systems (|Bode et all 2010l: Farris et al] 
20101: IO , Neill et al.ll2009l : IZanotti et all 120101 : iBode et al.l 1201 ll : 
Farris et al.ll201ll) . Postponing a more detailed calculation of this 
process to a future work, we have here shown that the calculation of 
the bremsstrahlung luminosities adopted in the above works leads 
to optimistic estimates, which should be regarded as upper limits. 
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APPENDIX A: EXTENDED GEOMETRIZED SYSTEM OF 
UNITS 

We recall that the definition of geometric units of time and lengths 
is obtained by setting the speed of light c and the gravitational con- 
stant G to pure numbers. This implies that seconds and grams of 
the cgs system can be written as 



Is 



2.997924xl0 10 ( - ) cm 



L0" 29 ( ^ ] 



(Al) 
(A2) 



Within this general setup, a convenient unit of space is still re- 
quired. The cm is of course a bad choice and the gravitational ra- 
dius r 9 = GM/c 2 is instead chosen. In order for this new unit to 
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be convenient with respect to the centimeter, the mass M of the 
system has to be sufficiently large. From the physical value of the 
solar mass and from (IA2t we find the relation between the cgs units 
and the new unit of length r g 



lcm 
Is 
lg 



( ^ I 



2.030281x10° 



5.027854x10 



/ V M 

2x 'Mr. 



-34 / C 



G 



(A3) 
(A4) 
(A5) 



It is also useful to write explicitly the conversion of rest-mass den- 
sity and luminosity between the two systems, namely 



p cgs = 6.1776X10 1 



3.6292 x 10 



G 

c 2 

59 / G 



M© 
M 



(A6) 
(A7) 



where /? cgs and p geo (as well as L cgs and L geo ) are the pure num- 
bers expressing the mass density (as well as the luminosity) in the 
cgs system and in the geometrized system, respectively. In the tra- 
ditional geometrized system c and G are set equal to unity. How- 
ever, for specific physical applications where very low mass densi- 
ties are encountered, the corresponding value of p geo may become 
prohibitively small. For this reason, it is convenient to assume a 
smaller value of G, such as G = 10 -10 . 

For convenience, we report the Eddington luminosity and 
the Thomson scattering opacity of electrons in geometrized units, 
namely 



^Edd = 3.4636 x 10" 



Xe = 3.628 x 10 22 Gp g 



Gj \MqJ 
M 



(A8) 
(A9) 



The extension of the geometrized system of units to the tem- 
perature can be obtained by setting to a pure number any physical 
constant containing the temperature. In this paper we have chosen 
to set mp/ks — 1, where m p is the mass of the proton, while ks 
is the Boltzmann constant. In this way the temperature is a dimen- 
sionless quantity and the transformation of the temperature from 
the dimensionless values to Kelvin is given by 



T K = l.C 



x 10 13 T„ 



(A10) 



In these extended geometrized units the radiation constant a ra d 
4a /c becomes 



1 ( i) (i£ 



(All) 
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